SUBROUTINE GRAMSCHMIDT(CPSI,NSTATES)
                                   !Gram-Schmidt orthonormalization
parameter (NX = 10001)
implicit real*8 (a,b,d-h,o-z) !variables starting with this letters are real by definition
implicit complex*16 (c) !variables starting with this letters are real by definition
double complex CPSI(NX, NSTATES), CAK(NX),CGRAM(NSTATES)

!common / constants / dx, cdt, ci, Ry, a0, am1, am2 !constants that will be used in several subroutines
common / constants / dx, dt, img, Ry, a0, am, am1, am2, ec, pi

!    | J > = | J > - SUM_i [ | i > * < i | J >

DO J = 2, NSTATES
 DO I = 1, J - 1
  CSUM = (0.E0,0.E0)
  DO IX = 1, NX
   CSUM = CSUM + DCONJG(CPSI(IX,I)) * CPSI(IX,J)
  ENDDO
  CGRAM(I) = -DX * CSUM
 ENDDO

 DO IX = 1, NX
  CAK(IX) = 0.E0
  DO K = 1, J-1
   CAK(IX) = CAK(IX) + CGRAM(K) * CPSI(IX,K)
  ENDDO
 ENDDO
 DO IX = 1, NX
  CPSI(IX,J) = CPSI(IX,J) + CAK(IX)
 ENDDO

ENDDO

DO J = 1, NSTATES !Normalization of all states
 SUM = 0.E0
 DO IX = 1, NX
  SUM = SUM + CONJG(CPSI(IX,J)) * CPSI(IX,J)
 ENDDO
 ANORJ = DX * SUM
 ACONS = 1.E0 / SQRT(ANORJ)
 DO IX = 1, NX
  CPSI(IX,J) = ACONS * CPSI(IX,J)
 ENDDO
ENDDO

return
end
